Fractional quantum oscillator and disorder in the vibrational spectra

We study the role of disorder in the vibration spectra of molecules and atoms in solids. This disorder may be described phenomenologically by a fractional generalization of ordinary quantum-mechanical oscillator problem. To be specific, this is accomplished by the introduction of a so-called fractional Laplacian (Riesz fractional derivative) to the Scrödinger equation with three-dimensional (3D) quadratic potential. To solve the obtained 3D spectral problem, we pass to the momentum space, where the problem simplifies greatly as fractional Laplacian becomes simply \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k^\mu $$\end{document}kμ, k is a modulus of the momentum vector and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document}μ is Lévy index, characterizing the degree of disorder. In this case, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu \rightarrow 0$$\end{document}μ→0 corresponds to the strongest disorder, while \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu \rightarrow 2$$\end{document}μ→2 to the weakest so that the case \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu =2$$\end{document}μ=2 corresponds to “ordinary” (i.e. that without fractional derivatives) 3D quantum harmonic oscillator. Combining analytical (variational) and numerical methods, we have shown that in the fractional (disordered) 3D oscillator problem, the famous orbital momentum degeneracy is lifted so that its energy starts to depend on orbital quantum number l. These features can have a strong impact on the physical properties of many solids, ranging from multiferroics to oxide heterostructures, which, in turn, are usable in modern microelectronic devices.

Although disorder (especially strong like substance amorphization) in a solid is usually considered to be a trouble for its possible device applications, its constructive properties have become increasingly clear in recent years as they give an additional possibility to fine tune their physical properties. Specifically, by varying the type and concentration of different imperfections (like point or extended defects), we can customize the characteristics of such material to meet specific requirements, needed, for instance, in electronics. Disordered semiconductors and dielectrics like organometallic halide perovskites 1,2 are widely used in photovoltaic cells, light-emitting devices, and nanolasers 3,4 . The above disorder influences phonon and electron spectra of a substance, leading to the distribution of the internal magnetic, electric and elastic fields. This is the case for so-called multiferroics, where ferroelectric and magnetic orders coexist 5 . It is well-known (see, e.g. 6,7 ) that strong disorder is well described by Non-Gaussian probability distributions, having so-called long tails. In other words, these distributions decay at infinities (sometimes much) slower, that Gaussian one. The most prominent example here is so-called Lévy flights [8][9][10][11][12][13] , which are Markovian random processes whose probability density functions (pdf) are Lévy stable laws, characterized by the Lévy index 0 < µ ≤ 2 . In the infinite space, it is profitable to define the corresponding pdf in terms of its characteristic function, i.e. Fourier transform. The pdf of a typical Lévy flight is usually determined by the fractional Fokker-Planck (FP) Eq. ( 8 ). In dimensionless units (diffusion coefficient and particle mass are set to unity) it reads where |�| µ/2 is an n-dimensional fractional Laplacian (Riesz fractional derivative 6,7 ) (1) ∂ t f (r, t) = −|�| µ/2 f (r, t) + ∇ ∇U(r)f (r, t)

The spectral problem for fractional quantum harmonic oscillator
In the units, where = m = ω = 1 (m and ω are the mass and frequency of oscillating particle respectively), the spectral problem for 3D fractional oscillator reads where the operator |�| µ/2 is defined by the expression (2) and ψ nlmµ (x) is the eigenfunction of a fractional quantum harmonic oscillator having the eigenenergy E nlµ for any specific µ value. In general 3D case, as it is typical for the quantum-mechanical problems in central force potentials, the eigenfunction ψ and eigenenergy E of the problem (4) depend on three quantum numbers n, l, m which are principal, orbital and magnetic quantum numbers respectively 18,19 . As we discussed above, the most profitable method to solve this problem is to pass to the momentum space, which permits to reduce the integral equation (4) to the differential one. Namely, as it is well-known that the combinations x 2 ψ(x, y, z) , y 2 ψ(x, y, z) and z 2 ψ(x, y, z) (here we suppress the indices nlmµ for brevity) in momentum space render as second derivatives ∂ 2 ∂k 2 x ψ(k x , k y , k z ) , ∂ 2 ∂k 2 y ψ(k x , k y , k z ) and ∂ 2 At the same time, the fractional Laplacian in momentum space becomes simply k µ , where k ≡= k 2 x + k 2 y + k 2 z is a vector k modulus. Latter definition immediately implies that 0 ≤ k < ∞ . This means that in momentum space, the equation (4) renders as where H k is the system Hamiltonian in momentum space. One note is in place here. Namely, the equation (5) shows that at µ < 2 , the famous orbital 3D oscillator degeneracy (see Ref. 18 for details) is lifted so that the eigenenergy E starts to depend on the orbital quantum number l. At the same time, as in the fractional case of µ < 2 , the time-reversal symmetry remains (it is lifted, for instance, in the external magnetic field 18 ), we still have degeneracy with respect to the quantum number m. Namely, the spectrum is invariant as m ↔ −m . That is why, in equation (5), we put index m explicitly in the wave functions, which are different for different m, although the energy E does not depend on m.
As for any central potential the wave functions are invariant with respect to rotation around the origin 18,19 , this is also the case for our 3D fractional oscillator. This implies that the wave functions ψ nlmµ (k) admit the separation of angular and radial (in k space) variables in the form where Y lm (θ, ϕ) = A lm P |m| l (cos θ)e imϕ are spherical harmonics 14,18 . Here P |m| l (cos θ) are associate Legendre polynomials of the first kind and A lm is the corresponding normalization coefficient, which is obtained from the following normalization condition 14,18 where d� = sin θdθdϕ.
It is well-known from "ordinary" quantum mechanics 18,19 that 3D quantum oscillator problem for µ = 2 admits the factorization of the wave function ψ(x, y, z) = ψ x (x)ψ y (y)ψ z (z) (here we also suppress the lower indices for a moment), which permits to reduce the problem to the three equations for 1D oscillator, containing www.nature.com/scientificreports/ three different energies, say E x , E y and E z . In this case, the initial energy E = E x + E y + E z . After problem being in such a way factorized, the eigenenergy can be derived in the well-known form 18,19 where n = n x + n y + n z with n x,y,z being the quantum numbers, corresponding to E x,y,z respectively. Such factorization in the "ordinary" case µ = 2 reduces the problem of the 3D oscillator to the combinatorial problem with respect to its 1D "components". In other words, the quantum numbers nlm of the resulting 3D problem comprise different combinations (obeying certain combinatorial rules, see Refs. 18,19 for details) of the initial n x , n y and n z numbers. It had been shown in 19 , that the possibility of such factorization at µ = 2 is based on additivity of both Laplacian operator (kinetic energy) and the potential ∼ (x 2 + y 2 + z 2 ) . At the same time, in r space (Eq.4), the fractional Laplacian (Eq. 2) is not additive. That is to say, it does not consist of a sum of the terms, containing separately x, y and z variables. This fact is "transferred" to k space, where the potential k µ ≡ (k 2 x + k 2 y + k 2 z ) µ/2 is factorable for µ = 2 (ordinary case) only. This signifies the difference between fractional µ < 2 and ordinary µ = 2 3D oscillator problems. This difference implies that for our fractional case we should proceed directly in 3D case. This, however, does not give an essential complication as problem still reduces to 1D one. Namely, after the decomposition (Eq. 6), we obtain following fractional 1D Scrödinger equation for the radial part R nlmµ (k) Here we also suppress the indices nlmµ in the functions R and energy E for a moment. In the quantum mechanics of the central force systems, it is also customary (see, e.g. 18,19 ) to pass to the function for which the Eq. (9) assumes particularly simple form which resembles that for 1D system. This form will be used below to construct trial functions for variational method as well as for numerical calculations, where the form (Eq. 11) permits to avoid divergencies at k = 0 . We note here, that if for the states with l = 0 , the function R(k) has a maximum at k = 0 , the corresponding function χ is zero by obvious reason.
As we mentioned above, the Eq. (11) can be readily solved numerically. On the other hand, even approximate analytical solution like variational one is much more desired as it can represent the useful orthogonal set, which can be further employed for the solutions of the problems, not obligatory dealing with fractional quantum mechanics. For instance, it could be well used to look for the solutions of the fractional FP equation (1) in the from of the expansion over the infinite sets of orthogonal functions. That is why we will solve the radial equation (11) variationally and confirm our results by direct numerical solution.
The Eqs. (9), (11) represent Schrödinger equation with the effective 1D potential |k| µ + l(l + 1)/k 2 . The plots of this potential at different µ and l are shown in Fig.1. It is seen from Fig.1a that for l = 0 the problem is identical  www.nature.com/scientificreports/ to that for 1D fractional oscillator 20 , while for l = 0 (Fig.1b), the centrifugal force drives it outwards. The only difference with the case of 1D oscillator is that now the problem is defined on positive semi-axis 0 < k < ∞ as k is now the modulus of a momentum vector. As it is seen from Fig.1b, the centrifugal term in the potential l(l + 1)/k 2 is dominant at k → 0 , while at k → ∞ the main contribution comes from k µ . This means that the whole spectrum of 3D fractional oscillator remains to be discrete similar to the "ordinary" case µ = 2 . Note that at small µ = 0.2 the potential begins to grow at very large k, while at µ = 0 the potential goes to 1 at k → ∞ . The same situation at l = 0 gives constant value one for µ = 0 , see Fig. 1a. Our numerical analysis shows that in 3D case despite the contribution of centrifugal term, at µ = 0 , the whole oscillator spectrum collapses into one value, see below. This property is similar to that in 1D case 20 . Note that at l = 0 and µ < 1 the potential at k → 0 the potential has nonanalytical behavior with infinite first derivative. At l = 0 this asymptotics is "killed" by the centrifugal term.

Variational solution. Classification of the states of 3D fractional oscillator
It is well-known that the variational principle of quantum mechanics 18 permits to find the approximate analytical expressions for the wave functions and eigenenergies of corresponding Schrödinger equations. This is especially true for the states of discrete spectrum, which are spatially localized. As usually, the "correct" (i.e. that which minimizes corresponding energy functional with minimal number of variational parameters) set of trial functions should be obtained on the class of functions, having the asymptotics, dictated by the corresponding differential equation.
Note that for the above variational principle to work, it is sufficient even to use the trial functions obeying the correct boundary conditions rather than have the "prescribed" asymptotics. In our case, this means that the trial functions have to be zero at k → ∞ . Such "simple" trial functions can be constructed on the base of the Gaussian one, inherent in an "ordinary" (i.e. that at µ = 2 ) 3D oscillator. In this case, the whole variational spectrum can be thought of in the form of "Hermite polynomials with unknown coefficients", where the latter coefficients are to be found from the orthogonality conditions, see below. In the next section, we shall compare the relative errors of the variational approaches, based on "fractional" (i.e. with predefined asymptotics) and "Gaussian" trial functions.
To construct the "fractional" trial functions, we should analyze the large k asymptotics of the χ(k) , given by the radial (in k space) Schrödinger equation (11). For that we observe that at large k we can neglect all the terms in the square brackets in Eq. (11) except k µ . This yields the equation for large k in the form d 2 χ(k)/dk 2 = k µ χ . The spatially decaying solution to this equation is proportional to Substitution of latter asymptotics into the above expression generates following wave function asymptotics at all admissible µ's Here we once more suppressed indices nlm. In the case of µ = 2 (ordinary quantum oscillator) the asymptotics (Eq. 12) reproduces well known result for "ordinary" 3D quantum oscillator 18 . It is seen that for all 0 < µ ≤ 2 the asymptotics (Eq. 12) decays sufficiently fast for corresponding integrals to be convergent. Namely, in worst case of k = 0 ψ(k → ∞) ∼ e −k √ 2 , i.e. decays exponentially. Good localization of the wave functions in k space yields their localization in coordinate space by Riemann-Lebesgue lemma, see e.g. 22 . Our analysis show that at small k the functions R(k) have the same asymptotics R ∼ k l (and accordingly χ ∼ k l+1 ) as in ordinary case, corresponding to µ = 2 18 .
To classify the states of the 3D fractional oscillator according to the quantum numbers nlm, we notice that although the Eq. (9) resembles that for the ordinary 3D oscillator, it is not equivalent to it as we have the term k µ instead of k 2 in the potential. The former term "spoils" the neat classification 18,19 of the states of 3D quantum mechanical oscillator so that in the fractional case µ < 2 there is no obvious way to introduce the principal quantum number n. We recollect that for ordinary 3D oscillator n = 2n r + l (see, e.g. 19 ), where n r is so-called radial quantum number, defining the number of nodes of the radial function for a given l. That is why to classify the states of the oscillator in our 3D fractional case we begin with general consideration of a quantum particle, moving in a central force field. Namely, to apply the oscillation theorem to the radial part χ(k) in Eq. (11), we arrange the eigenenergies corresponding to a given l in ascending order. It had been shown 18,19 , that the lowest eigenenergy for a given l corresponds to a nodeless function so that we assign n r = 0 to this state. The next (excited) state for the same l has one node, i.e. n r = 1 etc. This means that we may classify the states of our fractional oscillator by essentially two quantum numbers -n r and l. To be specific, for l = 0 (s -state; here we use the common notations of the states with the given l, see, e.g. 18 ) we have the whole sequence n r = 0, 1, 2, ... ; the same is true for l = 1, 2, 3, .. . Introducing the principal quantum number n = n r + 1 , we obtain following explicit classification of the states of the fractional quantum oscillator: 1s ( n r = 0 , hence n = 1 , l = 0 ), 1p ( n = 1 , l = 1 ), 1d ( n = 1 , l = 2),... 2s, 2p, 2d,... etc. In this case, the commonly known degeneracy of the states 1d and 2s, 1f and 2p etc, inherent in "ordinary" 3D oscillator, is lifted so that all states with different n r and l have different energies. This classification will be used below to construct the trial functions.
To construct the above functions, we need the asymptotics at small and large k as well as the information about the number of nodes, i.e. the numbers n r and l. Specifically, in terms of functions χ(k) , we have www.nature.com/scientificreports/ Here, A nl are normalization constants, found from the "1D-like" conditions Also, a nl are variational parameters, while the constants a i , b i and c i (here i = 1 − 4 ) are found from the orthogonality conditions where δ jj ′ is Cronecker delta. We note that the trial functions χ nl for different l (for instance χ 1p and χ 2d ) are orthogonal automatically because of orthogonality of spherical harmonics 14, 18,19 . That's why the integration in Eqs. (16) and (17) runs over k only. This implies, that only the functions with different n ′ s (and the same l) should be orthogonalized. Namely, in the particular case of the above trial functions we have explicitly The orthogonalization procedures Eqs. (17), (18) along with normalization of the functions χ nl determine uniquely the coefficients a i through variational parameters a nl . We have explicitly Here Ŵ(x) is the Ŵ-function 14 . The other normalization coefficients can be calculated, but their explicit expressions become progressively more cumbersome. Orthogonality condition < χ 1s χ 2s >= 0 yields (13) where is determined by the expression (20) and

Scientific Reports
Now, after variational determination of the parameter a 1s (see below), the parameter (20) starts to depend on the single parameter a 2s , which can also be determined variationally. Having determined the parameters a 1s and a 2s , we then use the orthogonality condition < χ 1s χ 3s >= 0 to express both c 2 and b 2 through a 2 . After substitution of the above found a 1s and a 2s into the corresponding normalization coefficient, the function χ 3s starts to depend on a 3s only. Such "chain rule" can be applied to determine uniquely all higher order functions, corresponding to l = 0 (s-state), l = 1 (p-state) and for any arbitrary n and l.
The expressions (13)-(15) show the algorithm of construction of trial functions. This algorithm is essentially the same as that for "ordinary" case of µ = 2 . Namely, the arbitrary trial function χ nl is the product of the factors k l+1 (small k asymptotics), e −a nl k 1+µ/2 (large k asymptotics) and a polynomial of degree n r = n − 1 with unknown coefficients. These coefficients, in turn, are determined from the orthogonality conditions (17) by the above "chain rule".
As usually, the variational solution of the Schrödinger equations should minimize the energies where H k is the Hamiltonian (5) and W nlµ are the variational approximations of the eigenenergies E nlµ . Normally W nlµ ≥ E nlµ . As the radial parts R nlmµ (k) of the wave functions (6) are well localized (see asymptotics (12)), the expression (23) could be rendered (by the integration by parts) to the more convenient form where R ′ = dR/dk . After substitution of the trial functions (13), (14) and (15) into the integral (24), with respect to the above "chain rule", we find the variational parameters a nl from its minimization. Such a procedure can be done for all wave functions of higher excited states, giving the approximate spectrum of the operator (5). Note that the substitution of found a 1s into W 1sµ (24) gives the approximate value of the ground state energy for all µ , the same with, say, a 2p gives the energy of the excited state W 2pµ and so on for higher eigenenergies. Although the expressions for higher excited states become extremely cumbersome (see above), the variational method permits to obtain the approximate analytical expressions for the eigenvalues and eigenfunctions of the operator (Eq. 5) for all admissible 0 < µ ≤ 2.
Substitution of the trial function χ 1s into the integral (Eq. 24) with subsequent minimization over a 1s yields which is similar to the expression for 1D fractional oscillator 20 . Further substitution of this value to the expression (24) generates the approximate value of the ground state energy for arbitrary µ It is seen that (W 1sµ ) min gives correct value 3/2 of the ground state energy for µ = 2 , corresponding to the ordinary quantum oscillator with the spectrum E n = n + 3/2 in our units. Below we shall see, that for the case µ = 0 all the spectrum shrinks into a single value E 0 = 1/2 , which is also obtained correctly from the expression (26). The dependence (Eq. 26) will be plotted below and compared with the numerical solution.
The same procedure with χ 1pµ gives that (a 1pµ ) min = (a 1sµ ) min , which is given by Ex. (25). The variational expression for the energy of the 1p state reads www.nature.com/scientificreports/ The expression (27) shows that at µ = 2 we obtain the correct answer 5/2. We recollect here that the energy of "ordinary" 3D quantum oscillator is determined by the expression (8), where n = 2n r + l , see above. As for 1p state we have n r = 0 (no nodes) and l = 1 so that n = 1 and the energy E 1p (µ = 2) = 1 + 3/2 = 5/2 . At the same time, at µ → 0 we have removable divergence lim µ→0 (2µ −µ/2 ) = lim µ→0 exp(−µ ln(2µ)) = 1 , while the rest of the expression (27) gives E 1p,µ=0 = 6Ŵ(3)/ Ŵ(5) = 1/2 , i.e. once more the correct answer. The dependence (27) will also be plotted below and compared with numerical results. We will see that the expressions (26) and (27) give very good approximate expressions for corresponding energies of a 3D fractional quantum oscillator for all admissible µ . Within the suggested variational approach, the energy of any desired level E nlµ can be evaluated analytically, although the calculations for higher excited states become very cumbersome, see above.

Numerical analysis
The easiest way to solve our problem numerically is to use the Eq. (11) for the function χ(k) , which has only one divergency (in the centrifugal term l(l + 1)/k 2 ) ) at k = 0 . The numerical solution of the Eq. (11) should be augmented by the boundary conditions χ(0) = χ(∞) = 0 . As usually for boundary value problems, our Eq. (11) can be solved numerically by its reduction to the eigenproblem for the corresponding finite-dimensional matrix. In this case, the energies E are the eigenvalues of the above matrix and the eigenvectors of the latter comprise the wave functions of the problem in momentum space. The transition to the usual coordinate space is accomplished by 3D inverse Fourier transformation. We note that the satisfactory accuracy of the numerical solution is achieved for typical matrix dimensions 10000 × 10000, which makes the task quite computer intensive.
As we can reproduce numerically the spectrum of the problem for arbitrary µ , we are now in a position to compare the results of variational and numerical treatments. We begin with comparison of different classes of trial functions (see above), namely aforementioned "fractional" (13)-(15) and "Gaussian", which we are going to construct below. As the consideration of a ground (1s in our case) state wave functions give the lower limit of the relative errors, here we consider just this function for "Gaussian" case. In terms of χ functions (10), we have where "G" in lower index stands for "Gaussian", A is normalization coefficient and σ is variational parameter. We note that at µ = 2 the function χ G1s is equivalent to χ 1s (see (13)) with a 1s = 1/(2σ 2 ).
Having both "Gaussian" (31) and "fractional" (26) and (27) expressions for variational energy levels of our 3D oscillator, we are now in a position to compare them with corresponding numerical values. Such comparison is reported in Fig. 2a,c. As is the case for the variational method, in both panels the variational curves lie higher than numerical ones as variational energy should be larger than its exact (in our case numerical) value 18 . It is seen that the agreement is a little worse for the 2s state, which has one node. For the nodeless states with n = 1 , the variational and numerical curves are indistinguishable in the scale of the plots. The quantitative error analysis for 1s state is reported in Fig. 2c. It shows that as both "Gaussian" and "Fractional" trial functions give exact results at µ = 0 and µ = 2 , the relative error should be maximal somewhere inside this interval. The inset in Fig. 2c shows that the maximal relative error occurs at µ ≈ 0.5 and does not exceed 0.6%. Latter error occurs for "Gaussian" trial function. This is because the latter function does not take into account the actual asymptotics (Eq. 12) of the wave functions.
The selected numerical wave functions in k-space are reported in Fig. 2b,d for different values of µ . It can be checked that functions are normalized, i.e. they obey condition (16). Note that corresponding "fractional" wave functions (13) are almost indistinguishable from numerical ones in the scale of the plot in Fig. 2b. At the same time, in the Fig. 2d, we choose the vertical scale so that the (minute) differences between "fractional" (variational)  Fig. 2 shows that functions for different µ have different decay rates, which are dictated by the asymptotics (Eq. 12). For instance, the 2s function for µ = 0.4 at k = 5 is not yet achieved its final part (compare to that for µ = 1.4 ), which occurs at k ≈ 15 . For smaller µ , the spatial extension of the wave functions in k space increases drastically. The functions for µ → 0 in k space are almost delocalized, which means that in coordinate space they are Dirac δ-function-like. This behavior is qualitatively similar to that of 1D fractional quantum oscillator 20 as well as for 2D 23,24 and 3D 25 fractional hydrogenic problems. The same tendency is seen in Fig. 2d. Here, the "Gaussian" trial function has stronger deviation from the numerical curve, while the "fractional" one goes closer to it. This difference is more visible at µ = 0.4 , which is close to the above value µ ≈ 0.5 , where the relative error in the energy is maximal. As numerical and variational curves are identical for µ = 2 (the corresponding trial functions are exact), the relative error is reported for µ = 1.4 and 0.4 only. The behavior of relative error reflects the well-known fact that different classes of trial functions can give a very good accuracy for the energy, while the approximation of the numerical wave function can be much worse. We note here that the values of relative error at both k → 0 and k → ∞ are not that reliable as in this case we are dealing with the difference of small (in the sense of machine precision) numbers, divided by those (small) numbers. That's why we limit ourselves by k ≈ 3 as for k > 3 because of latter effect, the error grows to infinity. It is seen from the insets in Fig. 2d, that in the intermediate region of k values, the maximal relative error of the "Gaussian" trial function is around 10%, while that for "fractional" one is less than 1%. Our analysis shows that for the "fractional" trial functions with nodes (including 2s and higher) the average error is 2-3%, the larger error around 3% occurs for the higher excited states with many nodes. This shows that the One more observation is in place here. Namely, the eigenenergies in fractional case µ < 2 are smaller than those for "ordinary" case µ = 2 . This means that in the system, where the substitution of conventional Laplacian by the fractional one in Schrödinger equation is admissible (like in strongly disordered solids like amorphous ones, see, e.g. 23 for discussion), it is energetically favorable for an oscillator to "become fractional".
As we have discussed above, the wave functions in the coordinate space can be obtained by the inverse Fourier transform. They are reported in the Fig. 3 for 1s (panel (a)) and 2s (panel (b)) states. It is seen from the Figure, that the curves are qualitatively similar to those in the momentum space. The regularities here are similar to those for 1D case 20 , i.e. the slowest decay have the functions for µ → 2 ("ordinary" quantum oscillator), while the fastest (like Dirac δ , see above) have the functions for µ → 0 ( µ = 0.4 in Fig. 3), corresponding physically to the strongest disorder. This may be one more signature of the famous Anderson localization 26 phenomenon.

Possible experimental applications. Outlook
The obtained spectrum of the fractional quantum 3D oscillator can be used for the calculation of arbitrary observable characteristic of real physical systems. As an illustration, here we consider the so-called oscillator strength or dipole matrix element (the matrix element of a dipole moment operator d = er ) for the transition between the energy levels with different orbital quantum numbers l and l ′ = l ± 1 . The latter relation stems from the selection rule, which gives zero (because of angular integration) for any combination of orbital quantum numbers, other than the latter one. The oscillator strength is important for the evaluation of such experimentally observable characteristics like the intensities of electric dipole transitions, which are the dominant effects of the interaction of a charge carrier in solids (electrons or holes) with the external electromagnetic field. One of the most important matrix elements of this type defines a transition between the states with l = 1 and l = 0 . As here we are working in momentum space, the formula for the oscillator strength contains the matrix element �a|k|b� for the transition between states a and b. In our dimensionless units, the expression for the oscillator strength reads (see, e.g. 19,27 ) where the transition occurs between states a (like 1p) and b (like 1s) with the eigenenergies E a and E b respectively. For the transition between the states 1p and 1s, the explicit expressions for matrix elements assume the form www.nature.com/scientificreports/ where Substitution of the matrix elements (Eqs. 33a-33c) with respect to numerically calculated R 00 and R 11 wave functions into the expression (32), generates the Lévy index dependence of oscillator strength components f i ≡ f i,1s−1p ( i = x, y, z ), which is reported in Fig. 4. It is seen that the oscillator strength goes to zero in the case of strongest disorder, corresponding to µ = 0 , see Ref. 23 for details. We observe (see above) that as µ → 0 , both E 1p and E 1s go to the same value (as, actually, the whole spectrum) 1/2. This means that for µ = 0 , the denominator of the expression (32) becomes zero. Despite that, the whole expression (32) tends to zero because of much faster decay of dipole matrix elements (33a)-(33c) in this case. The origin of this stems from coordinate space, where the wave functions of all states look very much like Dirac δ-functions. This means that the probability of optical transitions is much higher in relatively ordered regions, where µ > 1 . This fact may be useful in the theoretical description of the experimental data either in polymers, where chains can vibrate, or in disordered solids. Namely, the developed formalism can be well applied to the calculations of the properties of real physical systems, where disorder (like lattice imperfections and/or impurities) influences phonon and electron spectra of a substance, leading to non-Gaussian distribution of the internal electric, magnetic and elastic fields. A compelling example is the properties of the above multiferroics 5 . The non-Gaussian statistics due to disorder and frustration plays an important role in these substances [28][29][30] and it is challenging to apply the above formalism to explain their disorder (in phonon and magnon subsystems) phenomenologically in terms of introduction fractional derivatives to the corresponding Schrödinger equations. The dynamics of such "disordered phonons and magnons" can be described in terms of the above FP equation. Here, essentially the same variational approach can be used as there exist a mapping between FP and Schrödinger equation in terms of so-called Lévy-Schrödinger semigroup 31,32 . Of course, direct numerical solutions are always possible.
In this context, we also mention one more interesting physical problem regarding the onset of chaos 33 in the systems (like semiconductors 1,2 and trapped cold atomic gases 34 ) with the joint effect of the Lorentz force, Zeeman splitting, and spin-orbit coupling 35,36 . This suggests one more generalization of our 3D fractional quantum-mechanical oscillator problem. Namely, the spin-orbit interaction term can be added to the fractional Schrödinger equation (4). In this case, the solution will be more sophisticated as the wave function will be spinor now (see, e.g. 37 ) although the problem can be solved in the momentum space similar to the present case. This problem turns out to be very important for amorphous semiconductors with atomic vibrations 1,2,38 , where chaos can adversely influence possible optoelectronic, spintronic, and/or photovoltaic devices functionality.
Many physical problems can be reduced to that of a quantum harmonic oscillator. One of them is related to the theoretical studies of quark-antiquark bound states (so-called quarkonium, see, e.g. Ref. 39 ), which had been a hot topic several decades ago. In the context of 1D fractional quantum oscillator, the problem had been introduced by Laskin 17 , where the 1D potentials of the form |x| γ , 1 < γ ≤ 2 had been considered instead of quadratic one. At the same time, as the consideration was phenomenological by its nature, the linear potential had also been applied for that problem 39 . As it is well-known 18,19 , such potential in the 1D Schrödinger equation admits exact solution in terms of Airy functions 14 . In our context, such problem arises in Eq. (11) for l = 0 and µ = 1 . This problem had been considered by us earlier 40 in the context of Lévy flights confinement by linear potential. The  www.nature.com/scientificreports/ energies and wave functions of s-states (i.e. those for l = 0 ) are in a very good agreement with those obtained from exact solutions in terms of Airy functions. A quantitative comparison is out of the frames of the present work. Note also, that 1D and 2D problems for the fractional Schrödinger equation in a linear potential had been considered in the paper 41 . For 1D case this paper complements both our studies of 1D fractional oscillator 20 for the entire µ domain and those in the paper 40 for µ = 1 , where the exact solution in terms of Airy functions had been presented. The consideration of Cauchy distribution time evolution, made in 41 is also complementary to our study 42 . One more promising approach to non-relativistic quantum problems in central force potentials like r q , is to suggest an approximation to the spectrum of the problem for arbitrary q. This has been done in the paper 43 for the case of ordinary (i.e. non-fractional) laplacians in the Schrödinger equation. Moreover, the author 43 was instrumental to generalize the approximation both to the case of logarithmic potential, corresponding to (dr q /dq) q=0 = ln r and negative q like hydrogenic problems with q = −1 . It would be interesting in future to generalize the approach 43 to the case of fractional laplacians in corresponding Schrödinger equations.
Several one-dimensional problems (like fractional quantum harmonic oscillator and infinite quantum well as a limiting case of the potentials |x| p as p → ∞ ) had been considered in Ref. 44 . There, the corresponding timedependent Schrödinger equation had been solved by the matrix method. Different nonlocal operators, based on fractional Riesz derivatives had been discussed at length. Numerical and variational solutions had been presented. The variational solution of 1D fractional quantum oscillator problem had been reported in terms of Hermite and Laguerre polynomial type functions. This approach is complementary to ours in the sense that here we considered the "Gaussian" trial functions for comparison. The approach of the paper 44 is also complementary to our 1D case, considered in Ref. 20 . The consideration of the infinite potential well problem is complementary to our work 45 , where the same problem have been solved by the explicit matrix method. Note that the matrix method had been applied in our numerical solutions of 3D (present consideration) and 1D fractional oscillators problems.
In summary, we have studied the spectral problem for a 3D fractional quantum harmonic oscillator with arbitrary Lévy index 0 < µ ≤ 2 . Our main supposition here is that Laskin's construction of path integrals with Lévy measure is equivalent to "extraction" of probability density function from stochastic fractional Langevin equation and, in turn, to the assumptions made in seminal Anderson paper 26 . This permits us to assert that fractional Schrödinger equation accounts for disorder phenomenologically with Lévy index µ being the measure of the degree of disorder. The presence of potential in the system "tames" the initial Lévy distribution, making it decay faster than that in a corresponding free problem. In other words, here once more we have an interplay between the width of disorder distribution and system potential, which makes the square of the modulus of the corresponding wave function decay faster in space. This makes the problem of a fractional 3D oscillator resemblant to Anderson localization 26 in disordered systems.

Methods
The details of our theoretical methodology and those of working with fractional derivatives and fractional Laplacians, in particular, have been described in the sections "The spectral problem for fractional quantum harmonic oscillator" and "Variational solution. Classification of the states of 3D fractional oscillator". The numerical solutions of boundary value problems for partial differential equation (11) have been conducted using the commercial Mathematica software package.

Data availability
The datasets used and/or analysed during the current study available from the corresponding author on reasonable request.